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ABSTRACT 

The interaction between a massive binary and a non-self-gravitating circumbi- 
nary accretion disc is considered. The shape of the stationary twisted disc pro- 
duced by the binary is calculated. It is shown that the inner part of the disc 
must lie in the binary orbital plane for any value of the ivi@tfa.fy.chalmers.se 
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viscosity. 
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When the inner disc midplane is aligned with the binary orbital plane 
on the scales of interest and it rotates in the same sense as the binary, the 
modification to the disc structure and the rate of decay of the binary orbit, 
assumed circular, due to tidal exchange of angular momentum with the disc, 
are calculated. It is shown that the modified disc structure is well described 
by a self-similar solution of the non-linear diffusion equation governing the 
evolution of the disc surface density. The calculated time scale for decay of 
the binary orbit is always smaller than "accretion" time scale t acc — m/M (to 
is the mass of secondary component, and M is the disc accretion rate), and is 
determined by ratio of secondary mass to, assumed to be much smaller than 
the primary mass, the disc mass inside the initial binary orbit, and the form 
of viscosity in the disc. 

Key words: accretion, accretion discs-hydrodynamics-black hole physics- 
galactic centers-disc satellite interactions 
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1 INTRODUCTION 

In the present paper we discuss effects of interaction between a massive binary and circumbi- 
nary accretion disc. We consider in detail the evolution of the binary orbital parameters and 
the disc surface density. 

A massive binary surrounded by a gas disc may be found in different astrophysical sit- 
uations. For example, such a system may contain two pre-main-sequence stars surrounded 
by a protostellar disc. The tidal influence of this disc on the binary ormit may cause orbital 
migration, and in this way explain the observed small separations of some pre-main-sequence 
binaries. Similar processes may operate in protoplanetary discs producing short period plan- 
ets. The influence of the binary on the disc may be to significantly modify the disc structure 
by opening up a gap. The effect of such a gap may may be to reduce the total luminosity 
of the disc, and change its spectrum, and therefore it has observational consequences (see, 
e.g. Jensen et al 1996; Jensen & Mathieu 1997 for recent observations and an overview of 
theoretical models). 

Probably an even more important example is a supermassive binary black hole which is 
expected to be found in the centers of some galaxies. Such a system may be formed as a 
result of a galaxy merging process, and the evolution may end with the coalescence of the 
black holes. The burst of gravitational waves generated during such a coalescence should be 
easily detectable by the future generation of gravitational wave antennas. Just as in the case 
of binary stars, the interaction of a binary black hole with an accretion disc has prominent 
evolutionary aspects. The disc is able to collect and radiate away orbital binding energy, 
and thus takes part in the process of decay of the black hole separation distance. 

In fact the disc orbit interaction may determine the evolution of the binary separation 
distance over a definite range of scales. To see this let us briefly review the possible causes 
of evolution of the separation distance (Begelman, Blandford & Rees 1980, hereafter BBR). 
After the merger of the galaxies, the holes sink toward the center of of the stellar distri- 
bution (central stellar cluster) due to dynamical friction. When they are within a small 
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distance from the center ~ 1 — lOpc they form a gravitationally bound binary. After that 
the binary separation distance r& shrinks as a result of dynamical friction to another smaller 
characteristic scale r crit ~ — lpc where the binary becomes "hard" with respect to 
stellar population, i.e. the velocity obtained by a star during some typical encounter with 
the binary becomes larger than escape velocity of stars from the stellar cluster. Thus, a "loss 
cone" is formed in the velocity distribution of the stars in the cluster, and the binary can 
effectively interact only with stars which penetrate into "loss cone" at a rate determined by 
star-star relaxation processes. 

For separations ~ 10~ 2 — lpc decay of the binary orbit due to emission of gravitational 
waves appears to be insignificant, and the decay time scale is mainly determined by the 
two-body relaxation time of the central stellar cluster, which may be much larger than the 
cosmological time (BBR, see also Polnarev & Rees 1994; Quinlan 1996; Quinlan & Hernquist 
1997 and references therein for a more recent discussion). 

Note that the successive encounters of the stars with the binary lead also to a "random 
walk" of the binary center of mass around the center of stellar distribution (Bahcall & Wolf 
1976; Young 1977), and this effect may "smear out" the loss cone (Young 1977; Quinlan & 
Hernquist 1997). Although the problem of speeding up the orbital evolution of the binary 
might be solved in such a way, it would be very plausible to suppose there is another 
mechanism that can produce binary orbital energy loss at "intermediate" scales ~ 10~ 2 — lpc. 
In this connection we note that the size of the thin accretion disc around the more massive 
black hole (primary) may be as large as ~ 1CT 2 — 10 _1 pc, and therefore the interaction 
between the binary and the disc may "switch on" at the most critical "intermediate" scales 
~ 1(T 2 - KTV- 

The black hole binary-disc interaction can also have significant observational conse- 
quences in its own right. As we shall see in the next Section, during the initial stage of 
the binary-disc interaction the orbit of the smaller (secondary) black hole can be inclined 
with respect to the circumprimary disc, and in this case the secondary hits the disc twice 
during its orbital period. At the time of secondary-disc collision outflows of hot gas from 
the disc are formed, and the luminosity of the outflowing gas can be of order of secondary 
Eddington luminosity (Ivanov et al 1998, hereafter UN). Therefore the binary black hole 
may manifest itself as a periodic source of non-stationary radiation coming from the centers 
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of galaxies (see e.g. INN; Lehto & Valtonen 1996 and references therein). As we see below, 
during the final stage of the binary-disc interaction the binary orbital plane and the disc 
plane are aligned, and as for protostellar discs, the disc structure is modified by the gap 
formation process. The specific spectral features expected from such a modified disc can 
also be used to test the hypothesis of supermassive binary black holes. 

In this Paper we consider the evolutionary aspect of the interaction between a massive 
binary system (presumably, a binary black hole) and an accretion disc. We assume that 
our system consists of two companions with nonequivalent masses: a massive primary with 
mass M, and a secondary with mass m <C M; it also contains an accretion disc around 
the primary. We also assume that in course of previous evolution the binary separation 
distance r& becomes smaller than the disc size r^, so the disc mass inside the binary orbit 
Mdo(rb) < M^o(rrf), where Mdo(rd) is total mass of the disc (sometimes we shall call the 
distance the orbital distance of the secondary). We pay special attention to the case of a 
rather massive secondary: M^o (?"?>) < m < M^r^). Such a massive secondary significantly 
disturbs the inner part of the disc, and causes the evolution of the disc into a disturbed 
quasistationary state during a relatively short time-scale. The subsequent evolution of the 
binary parameters takes a much longer time, and can be calculated for a given perturbed 
state of the disc. 

In the beginning stages of the binary-disc interaction, the binary orbital plane is not 
necessary aligned with the disc plane, and therefore the gravitational field of the binary 
leads to disc distortions of two types: "horizontal" twisted distortions of the disc shape 
due to action of the non-spherical part of the binary gravitational potential, and "radial" 
distortions of the disc structure due to the tidal torque from the secondary acting on the 
disc gas particles which are orbiting close to the binary orbit. The "horizontal" distortions 
cause the formation of a quasi- stationary twisted disc where the inner part of the disc ( 
with radius smaller than some alignment radius r a i) lies in the binary orbital plane, and the 
orientation of the outer part of the disc is determined by the disc formation processes (see 
also Section 2,3). We calculate the shape of this twisted disc, and discuss the dependence 
of alignment scale on viscosity and mass ratio q = m/M (see also Kumar 1987,1990). The 
tidal torque of the secondary inhibits the disc gas from drifting inside the secondary orbit 
(Lin & Papaloizou 1979 ; Goldreich & Tremaine 1980), and so it accumulates in some region 
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outside the secondary orbit, and the radial structure of the disc in that region is modified 
(e.g. Syer & Clarke 1995). We calculate analytically the disc structure in that region in 
two ways: a) by using the laws of conservation of mass and angular momentum and b) 
by finding the appropriate self-similar solution of the diffusion equation which governs the 
viscous evolution of the surface density in the disc, and comparing the results with the 
results of simple numerical calculations. We also calculate the evolution of the secondary 
orbital separation r b (t) due to the torque acting between the modified part of the disc and 
the secondary. 

Our paper is organized as follows. In the next Section we estimate the characteristic time 
and spatial scales relevant for interaction between the binary and the disc, and qualitatively 
describe the possible evolution of such a system. All parameters are assumed to be typical 
for a supermassive binary black hole and galactic accretion disc. In Section 3 we calculate 
the shape of the stationary twisted disc around an inclined binary. In Section 4 we calculate 
the "radial" evolution of the binary and the disc. We summarize our results in Section 5. 
Although we mainly discuss the case of a supermassive binary black hole, our results can 
be applied to other possible situations such as stellar binaries or protoplanetary discs, with 
minor changes. The reader who is mainly interested in twisted discs around binary systems, 
can go directly to Section 3, and the reader who is mainly interested in the general aspects 
of the "radial" evolution of the binary and the disc, can jump to Section 4. 

We consider hereafter only circular binaries, and set eccentricity e to zero in our calcu- 
lations. 



2 MAIN PARAMETERS AND BASIC ESTIMATES 

In the beginning of this Section let us shortly review the basic properties of steady-state 
accretion discs. The detailed disc model is unimportant for us, and we employ simple a-disc 
models (Shakura & Sunyaev 1973) which are parameterized by the well known viscosity 
parameter a and accretion rate M. For a typical AGN disc the values of viscosity and 
dimensionless accretion rate are assumed to be small: a ~ 10~ 2 , M/Me ~ 10~ 2 , where Me 
is the Eddington accretion rate 

M B = ^~2xlO^M 8 (M / r ), 1 

CKt 
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here kt is Thompson opacity coefficient, and M 8 = M/10 8 M Q , and all other symbols have 
their usual meaning. We use later renormalized values of viscosity and accretion rate: a* = 
a/10" 2 , and M* = W 2 M/M E . 

All interesting disc quantities can easily be expressed in terms of disc semi-thickness h, 
and steady state surface density S . In the gas pressure dominated case, the semi-thickness 
is approximately proportional to the distance r, and to describe the vertical extent of the 
disc it is convenient to introduce the disc opening angle: 

6 = * « lO-V 171 ^" 1710 ^* 175 ^) 1720 , 2 
r 

where the distance r 3 = is assumed to be of order of the secondary orbital distance r b , 
and gravitational radius r g = 2GM/c 2 . [§] The steady state disc surface density is given by 
expression 

S « 10 5 «; 4 / 5 M 8 1/5 Mf 5 r- 3/5 g/cm 2 . 3 

With help of the expression (3) we can easily estimate the unperturbed disc mass Mdo inside 
the secondary orbit 

M d0 (r b ) « y7r£ (r 6 K 2 ~ 2 ■ lO^M^M^M^ 4 

The disc gas particles are orbiting around the central source with approximately Keplerian 
angular velocity 



Vk = \I^-* 0.7 • Mg-Vg 3/ V)~\ 



IGM 

and slowly drifting toward the center due to the action of viscous forces. The characteristic 
"viscous" time scale follows from eqs. (1,4): 

U = ^ ~ 10V 4/5 M 8 6/5 ilC 2 / 5 r 3 7 V 6 

Sometimes we also use the simple estimate for the drift time scale: 

The size of a thin 'standard' disc may be determined by the specific angular momentum 
of the gas which enters the disc, or by processes of star formation in the outermost regions 
of the disc, or any processes which can 'turn off' accretion. For the "standard" accretion 

§ In this Section we assume that opacity law in the disc is determined by Thompson opacity. The expressions for other 
reasonable opacity laws give similar results (see eqs. (44-47) for a general case). 
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disc with value of viscosity parameter a ~ 10 2 , the Toomre stability parameter Q can be 
rather small: 

Q « ^ « ™ « 0.5«I/ 10 M 8 - 13 / 10 M- 2 / 5 r3 27/20 

In the disc region where Q < 1 the disc is gravitationally unstable and locally gas may 
collapse to form low-mass stars (e. g. Spitzer & Saslaw 1966, Illarionov & Romanova 1988). 
We do not consider such a disc here as we adopt disc parameters such that Q > 1. As 
mentioned above, we consider the case when the interaction between the secondary and the 
disc occurs on scales smaller than the disc size r^. That condition can be described in terms 
of characteristic radial scale r m -the scale where the disc mass is equal to the mass of the 
secondary: 

r m « 3 ■ 10 3 at/ 7 mfM 8 U/7 M^ 7 r 9 

where tuq = m/lO 6 M . The secondary interacts only with gas contained inside the scale 
r m , and the disc region outside this radius is in the unperturbed steady state 0j- The gas 



m . 

" ,:c ~~ M ' 



contained inside the radius r m is accreted during the time t 

tacc = t u0 (r m ) = 5 ■ 10 8 ^M~ V, 8 

Ms 

and this time is larger than any characteristic time scale relevant to our problem. In the 
case of a supermassive binary black hole, the 'accretion' time (8) must be smaller than 
cosmological time in order to provide the source of gravitational radiation (BBR). 

Now let us describe qualitatively the scenario for binary evolution due to the interaction 
with the accretion disc, and estimate the relevant characteristic spatial and time scales. 
As we will see the interaction of the secondary with the disc gas leads to alignment of the 
secondary orbital plane and the plane of the inner disc. After the disc alignment has occurred, 
the orbital separation starts to shrink in the radial direction. Therefore the evolution of the 
binary and disc can be roughly divided into two successive stages -the alignment stage and 
the stage of 'radial' evolution, and we consider these stages in turn. 

First we assume that the secondary orbital plane is inclined with respect to the disc 
plane. The effects of the interaction between the inclined secondary and the disc can be 

^ Actually, in case of a very massive secondary, the disc structure is close to the steady state solution at scales smaller than 
r m , see eqs. (63,64) below. 
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considered in two regimes. The first is when highly non-stationary distortions of the disc 
result from direct collisions between the secondary and disc (s.-d. collisions). The second is 
when the long-ranged, averaged over many orbital periods, interaction between the binary 
and the disc produces a twisted warped disc. 

As we will see the evolution of the secondary orbit due to s.-d. collisions is important 
for relatively low-mass secondaries, and long-ranged interaction is important in the opposite 
case. 

Let us estimate effects of the direct collisions between the secondary and the disc. The 
inclined secondary hits the disc twice during an orbital period £ or &, which is approximately 
equal to the Keplerian period: 

t or j, ~ 27rfi^ 1 . 9 

During a collision, the secondary significantly perturbs the disc gas at distances from it 

comparable with the "accretion" radius: 

2Gm 14 

3 x 10 m 6 r 3 cm, 10 



•a i — * — * i o 
\Vsec -Vd\ 2 

where v sec is the velocity of the secondary which intersects the disc at the "collision" radius 
r c ~ r b , v d is the disc velocity at the point of collision, and \v sec \ ~ \vd\ ~ r c VL{r c ). The 
perturbed gas is heated up by the shock induced in the disc. An amount of gas of mass 
~ vrSr^ is accreted by the secondary, another amount of approximately the same mass 
attains velocities greater than the escape velocity in the binary potential and leaves the 
system (UN). As result of the collision, the linear momentum of the secondary is changed, 
with the momentum change Sp sec being proportional to the total mass of perturbed region 
> 27rE(rb)r^, and it is directed opposite to the relative velocity of the secondary v sec — 
Therefore the change of velocity of the secondary Sv sec can be written as: 

0V S ec ~ {Vsec ~ V d ) , 11 

m 

here the correction factor A > 1 is mainly determined by the disc distortions at scales larger 
than r a The secondary orbital parameters change on account of the change of velocity, 
and this leads to the dragging of the secondary orbit into the disc plane. To see this, consider 

II Wc do not take into account the resonant excitation of waves in the disc due to the secondary passing through the disc. 
Estimates of this effect have been made e. g. by Artymowicz 1994. 
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the evolution of the inclination angle (3\ the angle between the specific angular momentum 
of the secondary l sec , and the specific angular momentum vector of the disc gas l^ = f c x v d 
(see Fig. 1 below). From eq. (11) it is easy to obtain the equation for the change of this 
angle during the collision |_1: 

6(3 = _ 2nZrlA \l d \ ^ 

m | hec | 

One can easily see from this expression that there is only one stable value of (3: (3 = 0. Thus, 
collisions between the secondary and the disc tend to align the angular momentum vectors of 
the disc and the secondary. When these collisions determine the evolution of the secondary 
orbital parameters, the secondary is progressively dragged into the disc plane until it is in a 
coplanar orbit. The associated timescale tdragg can be easily estimated with help of eq. (11) 
to be 

m 

tdrag ~ 2^rlA torb - 

As we mentioned above the b.h-disc collisions result in an outflow of disc gas at a rate 

M out > E r 2 a n K , 14 

Obviously, this outflow must be compensated by a radial inflow of mass through the disc, 
and therefore the surface density in eqs. (11-13) is close to the steady state value (3) only if 
the condition M out < M is satisfied. This condition gives: 

q < ^a^S ~ 1.5 • lO^Ba'J 2 , 15 

where black holes mass ratio 

m 

q = w 

the parameter B = 5/1CT 3 ~ 1, and we use the standard equation of stationary disc ac- 
cretion: M = 37rz/£o and kinematic viscosity v = a5 2 r 2 QK- In the opposite case q > aV 2 5, 
the secondary depletes gas from the region of the disc near its orbit, so the surface density 
in this region should be much less than the unperturbed surface density S . Now an upper 
limit for the surface density E(r c ) may be estimated by equating the rate of outflow of disc 

** In the analogous case of star-disc collision, similar calculations and references to earlier work can be found in a paper by 
Vokrouhlicky & Karas 1998. 
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gas ~ S(r c )r^fi^, and the disc accretion rate M, and that gives 



of 
q 2 



nr p ) ~ —S < S 16 



In the limit <C 1 the drag timescale is of order of 'accretion' time scale t acc . 

In the case of a sufficiently massive binary, the disc midplane is expected to be aligned 
with the binary orbital plane due to the influence of quadrupole component of the binary 
gravitational field. Schematically this process can be divided again into two successive stages. 
The planar disc first evolves into a quasi-stationary twisted configuration (the twisted disc) 
during some characteristic alignment time t a i. In the twisted disc, the inner part of the disc 
is aligned with the binary orbital plane, and the outer part does not change its orientation. 
The radius out to which the disc is aligned (r a j), as well as the alignment time t a i, are 
estimated in the next Section (see also Kumar 1987, 1990). Both quantities strongly depend 
on the value of viscosity parameter a. For sufficiently large values of a > 5, we obtain: 

r a ii ~ (aq) 1/2 5~ l r b , 17 

and diffusion-like decay to a quasi-stationary twisted disc proceeds during the time 

t a i ~ a\ ~ a5~ 2 Vf^. 18 

For small values of a, the alignment scale does not depend on viscosity (Ivanov & Illarionov 
1997, hereafter II) 

r a n ~ q 1/2 S^ 2 r b , 19 

and the evolution to the stationary configuration has wave-like character (Papaloizou & Lin 
1995, hereafter PL). Now the alignment time t a i is approximately equal to the 'sound' time 
t s : the time taken by sound wave to cross the radius r: 

t s ~ r^OaO- 20 

Obviously, the alignment scale must be greater than the binary orbital distance: max(r a n, 7^2) > 
r b . That condition gives: 

5 2 

q > — when a > 5; q > 5 when a < 5. 21 
a 

When the inequalities (21) are broken, the alignment effect is absent. 

After the formation of a quasistationary twisted disc, the orientation of the binary orbital 
plane itself slowly changes with time due to the influence of the gas accreting through the 
twisted disc. When the gas is accreting through the twisted disc, the orientation of its angular 
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momentum is changing, and therefore the orientation of the binary orbital plane must also 
be change with time due to the law of conservation of angular momentum (e. g. Rees 1978). 
This fact allows us to estimate the change with time of the binary inclination angle (3 out as: 

Pout ~ -^pi, p out (t) oc e~~ , 22 

where index out means that the inclination is determined with respect to the orientation 

of the disc at radii far from the binary. The characteristic time-scale of this effect £* is 

determined by the ratio of specific angular momentum of the binary to the specific angular 

momentum of the disc at radius r ~ r a i, and by the 'accretion' time t acc : 

. r b v d (r b ) 
r a iVd{rai) 

Here the alignment scale is given by eq. (17) and a is assumed to be rather large (a ~ 1) 
(the case of small a C 1 is more complicated and must be considered separately, see e.g. 
Scheuer & Feiler 1996). When the secondary counter rotates with respect to the disc gas ( 
retrograde rotation), we have ^ < 0, and the time t ev is formally negative. This means that 
the inclination of the binary orbital plane with respect to the disc plane grows with time, 
and the binary may turn over and become prograde during the time ~ t ev . Note, however 
that a highly inclined twisted disc may be unstable with respect to instabilities of a different 
kind, and so the evolution of a retrograde inclined binary is unclear. We will assume below 
that the binary is rotating prograde with respect to the disc gas. 

From eq. (23) and eqs. (18, 20) it follows that t ev ^> t a i when m 3> This inequality 
as well as the inequalities (21) tell us that the disc alignment effects are important only for 
sufficiently massive secondaries (typically q > 1CT 3 ). In the opposite case the secondary is 
quickly dragged into the disc midplane by the competing effects of secondary-disc collisions. 

Now let us assume that the binary and the disc are aligned at scales of interest, and 
consider the next 'radial' stage of the evolution of the binary and the disc. At first let us 
look at the case of the small mass of the secondary m < M<j. As we mentioned above a low 
massive secondary must lie in the disc plane corotating with the disc gas. Orbiting in the disc 
the secondary transfers the angular momentum to the neighboring gas particles which are 
orbiting upstream, and gain the angular momentum from downstream gas particles due to 
resonance interactions/ non-resonance scattering of these particles (Lin & Papaloizou 1979 
Goldreich & Treimain 1980). This can lead to the formation of gap around the secondary 

© RAS, MNRAS 000, 



12 P.B.Ivanov, J.C.B. Papaloizou and A. G. Polnarev 



which is free from the gas, provided the gap formation criterion is satisfied (see, e.g. Lin & 
Papaloizou 1979, 1986) 

q> wrK 24 

where is the binary separation distance. For a-discs the criterion (24) can be rewritten as 

q > ^-a5 2 ~ 3 x l(T 7 a*. 
8 

One can see from this relation that the criterion is obviously satisfied for any reasonable 
mass of the secondary. Lin and Papaloizou showed that when the gap is opened the orbit of 
secondary is shrinking during viscous time scale t u $ [n]. 

At the end of this section let us briefly discuss the second possibility m > ( this case 
is discussed in detail in Section 3). Now the radial drift of the secondary proceeds during 
the time much larger than t U Q. The secondary torque acting on the downstream gas particle 
prevents them to drift inside the secondary orbit, and the downstream gas is accumulated 
in some region outside the secondary orbit. The structure of the disc in that region is 
significantly modified. The radius out to which the disc is modified grows with time (see 
Section 3). The accumulated mass of the gas near the secondary orbit increases the torque 
acting on the secondary, and pushes the secondary toward the primary. The secondary orbit 
is shrinking during the characteristic evolution time-scale t ev (t u0 <C t ev < t acc ), which is 
evaluated in the Section 3. 



3 THE TWISTED ACCRETION DISC AROUND THE INCLINED 
BINARY 

Now let us consider in more detail the evolution and the structure of the twisted disc around 
the binary. At first we need to choose a convenient coordinate system. As a basic coordinate 
system we use the Cartesian coordinate system (x, y, z), which has its origin (O) at the center 
of mass of the binary. The plane (XOY) coincides with the binary orbital plane and the axis 
(OZ) is directed perpendicular to that plane. As it follows from our previous discussion the 
binary orbital plane is slowly evolving with time, but when considering the dynamics and 

tt Note that the secondary can shrink its orbit even though it is not massive enough to open a gap (Ward 1997 and references 
therein) 
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Figure 1. The r, ip, £ coordinate system shown with respect to the Cartesian coordinate system. P is the projection of the 
point R on to the disc ring plane. 

the structure of the twisted disc itself the effect of this evolution is unimportant. Since 

the twisted disc has some inclination with respect to the binary plane, and this inclination is 

different for different radii r, when studying the local properties of the disc it is convenient 

to use another locally cylindrical (twisted) coordinate system (r, ■?/>,£), (Petterson 1977). 

This system consists of rings centered at the point (O), which have different (with radius) 

inclination angles /3(r,t) with respect to the binary orbital plane, and the lines of nodes of 

these rings are rotated through some angle 7(7", t) with respect to the (OX)-axis, see Fig. 

1. The position of the gas elements on the ring is characterized by the angle ip, and the 

coordinate £ characterizes the vertical displacement of the gas element out of the ring plane 

(see Fig. 1). Assuming that the twisted disc is evolving on a time scale much larger than the 

orbital time scale, and the spatial scales of the disc twist are much larger than the binary 

orbital distance, we hereafter use the time-averaged over many orbital periods monopole and 

quadrupole components of the binary potential 0(r, ■?/>,£): 

GM Gmrt 3R 2 2^sin^ 

where the spherical coordinate R = y/r' 2 + £ 2 , and we assume (3 < 1 f% 

When disc twist is small < 1), trajectories of the gas particles in the twisted disc are 

« In the expression (25) the assumption m -C M can be relaxed. The potential (25) describes the time-averaged gravitational 
field of the circular binary with components having the comparable masses mi, m,2, provided M = mi +rri2, and m = J ^ 1 _^ 2 ■ 

13 
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ellipses with small eccentricity e. The influence of the non-spherical part of the potential (25) 
on these trajectories can be described as the following. Neglecting the interaction between 
gas particles and then considering the free motion of these particles in the potential (25), 
one can see that the quadrupole contribution to the potential causes the precession of the 
major axis of an elliptical orbit with the frequency: 

^ap = 7?(— ) 26 
4 r 

and precession of the lines of nodes with the frequency: 

Q = —Q 27 

The apsidal precession is prograde and the nodal precession is retrograde with respect to 
the orbital motion. As we will see, the signs of these precession frequencies determines the 
shape of the stationary twisted disc in the low- viscosity limit. 

The dynamics and the structure of a twisted disc can be described in terms of two 
complex variables: W(r, t) = /3e n , and A(r,t). The variable W determines the inclination 
and the rotation of the disc rings and the variable A(r, t) determines the deviation of the 
gas trajectories from the circle of the radius r. Namely, the eccentricity e and the position of 
the apsidal line of the elliptical trajectory ip can be expressed in terms of A: e = — 2-0#|A|, 
ipQ = arccos (Demianski & Ivanov 1997). In the low- viscosity limit a < 1, the dynamical 
twist equations can be written as 

A = - j^W jr + {iVL ap - an K )A, 28 

W = 7r(^ 2 A) + in np W. 29 

r or 

The derivation of these equations lies beyond the scope of our paper, but can be found 
in papers by Papaloizou and Pringle 1983, Papaloizou and Lin 1995, Ivanov and Illarionov 
1997, Demianski and Ivanov 1997. Note that the equations (28-29) are applicable not only 
to the case of a binary system, but also for general perturbing forces which can cause apsidal 
and nodal precession of the disc gas trajectories in the linear approximation over the angle 
(3 (II, Demianski & Ivanov 1997). 

We do not solve the time dependent form of these equations in our work, and restrict 

§§ The variable W is connected to the variable g used by PL by W = (ig)* ■ 
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ourselves to qualitative analysis only fi\ For simplicity in our analysis of the time dependent 
case we also neglect the terms proportional to the precession frequencies -the last terms on 
the right hand side of eqs. (28), (29) (for more detailed analysis see, e.g Papaloizou & Lin 
1997). These terms determine the shape of the stationary solutions of the twist equations, 
but are not important for the character of time behavior. Let us also look for the solution of 
eqs. (28), (29) in the WKB-limit: W oc ^M*)*-*^ kr < 1. We obtain the dispersion relation 
u(k) in the form (PL): 

1 



u = ~ (iaQ K ± yj k 2 v 2 - a 2 Q 2 K ) , 30 

where v s = 5rVtx is the sound velocity. In the case k > ^ eq. (30) describes the propagation 
of the waves with velocity vw ~ v s /2. In the opposite case k < the mode corresponding 
to the sign (+) in the eq. (28) describes the decay of the ellipticity of the gas trajectories 
due to the action of viscous forces. The mode corresponding to the sign (— ) is described by 
the relation 

8 2 

w « i — ( r 2 k 2 )n K . 31 
4a 

This relation implies that the relaxation to the stationary configuration has diffusion-like 
character, and the characteristic time scale of this effect t a i(k) ~ a<5 _2 (rfc) _2 t or ,{,. Considering 
the modes with k ~ r, it is easy to see that the transition from wave-like to diffusion-like 
propagation of the twist corresponds to a > 5, and when a > 5 the disc alignment time-scale 
t a i can be estimated as: 

ted ~ t a i(k ~ r _1 ) ~ a5~ 2 t orb . 

In the opposite case a < 5, we have: 

t a i ~ r/v s ~ 5~H orb . 

In stationary case eqs. (28, 29) can be written as: 

d *3/ 2 d _ w _ 3i| w ^ 



dx A dx x 5 / 2 
A = a — -^qx , 

where we use the explicit form for the precession frequencies (28, 29), introduce the dimen- 

^ For the discussion of numerical solutions of these equations see Kumar 1990, PL. 



© RAS, MNRAS 000, 



16 P.B.Ivanov, J.C.B. Papaloizou and A. G. Polnarev 



sionless distance x = -jr, and new parameters a = j and q = |. The eq. (32) is similar 
to the equation describing a stationary twisted disc around a Kerr black hole (II) except 
for the (— ) sign in front of the forcing term (the last term in (32)). Ivanov and Illarionov 
found radial oscillations of the inclination angle of a low-viscosity stationary twisted disc 
around a Kerr black hole, and suggested that this effect holds for any nearly Keplerian forces 
which cause apsidal and nodal precession of the disc trajectories provided that these preces- 
sions are in the same direction. (jf^ > 0). They also suggested that when the precessions 
counter-rotated < 0) the disc inclination must join smoothly to the orbital plane of the 
object. A twisted disc around a binary system gives a natural example of counter-rotating 
precessions, and we see below in our case that the disc lies in the orbital plane for any value 
of viscosity. 

Equation (32) cannot be solved analytically; its numerical solution is presented at Figs. 

2, 3. However, we can obtain approximate analytic solutions of this equation for the case 

a ^> 1, and a < 1. When d > 1 we can set A « a, and the approximate solution of eq. 

(32) is expressed in terms of Kelvin functions ker(x), kei(x): 

2 3/A W out i 

W = —j ou y£ (kerUyx) + ikeii{y x )), 33 
r (l) 

where T(x) is the gamma-function, W out = W(r — > oo), y\ = raU -, and 

r all = (3qa) 1/2 r b . 34 
Hereafter we assume W(r — > 0) —>■ as an inner boundary condition. When r — > 0, 

, , r all i ( r al\ i 3-7T nl \ 

W w C 1 (r)e~v^e~ l[ v^ + ^~ lout \ 35 

where C\(r) is a slowly changing real function. When a <C 1, A ~ — jiqx~ 2 , and the solution 
is expressed in terms of modified Bessel function K(x): 



where y% = (— ) 2 , and 



When r — > 0, 



(3g) 1/2 

r a i2 = — - — r b . 37 



W oc e-PP) . 38 
Comparing eqs (35), (38) we see that in the case of a relatively high viscosity, the disc rings 
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Figure 2. The dependence of the inclination angle /3(r)/f3 ou t on radius r/r;, for q = 100, and for the different values of 
a. Curves 1, 2 and 3 correspond to a = 0.01,1,100 respectively. We also plot two simple analytic approximations for these 

curves which correspond to small and large values of a. The dotted curve /3 = e~ ' r ) (r a (2 = 8 . 7r^ ) , and the dashed curve 

P = e ^ r (r , a (2( Cf = 100) = 173.2rj,) are the approximate solutions of the twist equation in the limit of high viscosity (see 
eqs. 35, 38). 
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Figure 3. The result of the integration of the twist equation (32) in a parametric form (\Pi = i?eW, *2 = IrnW). The larger 
values of /3 = yj Nl^ + correspond to the larger values of radial distance. The curve 1 correspond to the curve 1 of Fig. 2 
(q = 100 hereafter and a = 0.01), the curves 2 and 3 correspond to the curves 2,3 of Fig. 2 (a = 1 and a = 100 respectively). 
When viscosity is small, the rotational angle 7 = arctan (^) is not changed with the change of r, for larger values of viscosity 
the rotational angle is changing with the distance. 

are twisting with r, but in the low viscosity limit the effect is absent. This might give a 
possibility of testing the value of viscosity in circumbinary twisted discs. 

The general case a ~ 1 must be analyzed numerically. We plot the results of numerical 
integration of the twist equation (32) in Figs. 2,3. 
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As follows from Fig. 2., the alignment scale decreases with decrease of the viscosity 
parameter in agreement with our analytic estimates, and the expressions (33), (36) give a 
good approximation to the numerical solutions of eq. (30). Fig. 3 illustrates the change of 
the twist angle 7 with the distance r. One can see from this Figure that the angle 7 does 
not change significantly for small values of the viscosity parameter. 

In deriving eq. (32) we have neglected the self-gravity of the disc. The effect of self- 
gravity is manifest in two different ways: it gives an additional non-local forcing term F sg 
in eq. (32), as well as a contribution to the resonant denominator A. The term F sg must be 
compared with the term induced by pressure - F p - the first term on left-hand side of eq. 
(32). According to PL, the forcing term F sg can be estimated as: F sg ~ ~§^Fpi an d we can 
neglect this term assuming that viscosity takes a rather small value |jnt : 



On the other hand, the correction to the resonant denominator is more important in case of 
low viscosity, and can be estimated formally by setting a = 0. As well as giving a correction 
to the Newtonian potential due to the secondary, self-gravity induces apsidal precession of 
elliptical trajectories with a frequency il sg ~ ^VLk- This frequency must be compared to 
the frequency Q ap at characteristic scale r a i 2 : 



Thus, for estimating precession frequencies, the correction due to self-gravity appears to be 
less important than the correction due to gravitational field of the secondary if Q > 1. 

At the end of this Section we note that the low viscosity a < 1 twisted disc may be 
unstable with respect to shear instabilities if the "initial" inclination angle (3 out is sufficiently 
large (3 out > 5 (Coleman & Kumar 1993, II). It was suggested (II) that the shear instability 
might result in development of turbulence and so effectively increase the value of the viscosity 
parameter. In this case the effective scale of the disc alignment increases and equation (37) 
gives only a lower limit to the alignment scale. 



In a similar context a self-gravitating twisted disc with (Q ~ 1) has been recently discussed by Papaloizou, Terquem and 



a < 5Q. 



M d (r al2 ) 
M 




< Qap{r a i 2 ) ~ Sfl K . 



Lin 1998. 
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4 BINARY ORBIT EVOLUTION IN THE PLANE OF THE ALIGNED 
DISC 

Consider now in detail the radial evolution of the disc and the secondary. In this Section we 
assume that the disc plane is aligned with the orbital plane of the binary near the secondary 
orbit, and the secondary rotates in the same sense as the disc. Hereafter we adopt the 
following assumptions: a) we assume that the mass of the secondary is much smaller than 
the primary mass M, but larger than the characteristic disc mass b) the evolution due 
to interaction of the secondary and disc is switched on at some moment of time t — to when 
the secondary lies in the disc at some radius r ~ r c < r<j. The disc is taken to have a steady 
state structure at this moment c) Assuming that the inner (with respect to the secondary 
orbit) part of the disc is relatively quickly swallowed up by the primary, we do not consider 
the evolution of that part, d) We assume the coefficient of kinematic viscosity v to be a 
power-law function of cylindrical coordinate r and surface density S: 

v = kY, a r h . 39 

The expression (39) can be specified for a given model of the accretion disc. In the a-disc 
model the relation between surface density and viscosity follows from consideration of the 
energy balance. Assuming the power-law dependence of the opacity k on the density and 
temperature 

K = K 0P ^T^ 1 

a simple calculation gives (Lyubarskij & Shakura 1987) 

k = f{a 8+ ^' 2 ^4n 2(4 '^\GMY 2 ^ + ^^ 2 (ca r )- 2 y, 40 

and 

2(2 + /i X ) 3(4-/x 1 -2 M2 ) 

a = , o = , 

e 2e 

and e = 6 — 2/i2 + H\- Here a r is the radiation density constant, TZ is the gas constant, c 
is the speed of light, and the standard system of units is assumed. The parameter / ~ 1 
depends on the disc structure in vertical direction and is unimportant for us. If Thompson 
opacity dominates, k = ktk = 0A(cm 2 g~ l ), a = 2/3, and b = 1. If the opacity is determined 
by free-free processes, Kff ~ 0.8 ■ 10 23 pT~ 7 ^ 2 (cm 2 g _1 ), and then a = 3/7, b = 15/14. 

An analogous problem was formulated and discussed by Syer and Clarke 1995 (hereafter 
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SC). However, our results differ from those obtained by SC. As we see below the effect of the 
interaction between the secondary and the disc leads to accumulation of the disc gas in the 
region of the disc confined between the binary orbital radius r^, and some 'binary influence' 
radius r in fi(t) > r&. The radius r in fi grows with time. Syer and Clarke approximated r in fi as 
a constant, and this led to an overestimation of the disc surface density near the secondary 
orbit, and underestimation of the binary evolution time-scale in their calculations. 

4.1 Evolution of the disc and the binary 

The master equation for the evolution of the disc surface density X can be written as (e.g 
Lightman Sz Eardley 1974, Lyubarskii & Shakura 1987): 

9^ = z±d ( 1/2 d (rl/ ^ +a) ) _ 41 



dt r dr \ dr 

Strictly speaking this equation is not valid near the secondary orbit, where the tidal inter- 
action exchanges the angular momentum between the disc and the secondary. We assume 
hereafter that the tidal torque is concentrated in the narrow disc ring near the secondary 
orbit, and eq. (41) is valid at all distances outside this ring. Following this assumption the 
influence of the secondary on the disc can be approximated as an effective boundary condi- 
tion for eq. (41), which takes into account the exchange of the ^-component of the angular 

momentum between the secondary and the disc (SC) 

3 1 dr^ 

{27rr 2 £(-i/f2^- + Vdrrflx) H — mrVtx—r- } r =r b — 0, 42 



where all quantities are taken at r = r;,(i), and the secondary radial velocity -Q is of order 

of the disc gas drift velocity Vdr near the secondary orbit. In general, 

3 d 



Vdr 



(r 1/2 z/S) , 43 



YjT 1 ! 2 dr 

In the case of a massive secondary (m ^> Md), the second term in the round brackets can 
be neglected. 

Prior to the moment t = to, we set |j| = into eq. (41). In this steady state case we 
have 

M 

7]o = Z/ £ = 7^, 44 
37T 

1 

i M 1+a b 

S = (3Jfc)-3TS(_) r -—«, 45 



7T 
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1, , 1 ,M, 1 + a b-a-l 



V rO = -o( 3fc ) 1+a ( — ) r 1+ ° , 46 
Z 7T 



i ,M, 1+a 



Uo = -, r = 2(3fc) _ ^(— ) r~ % 47 

\V r0 \ 7T 

where 

c = 2(a + 1) -b. 

We assume below that c > 0. 

The presence of the secondary in the disc significantly influences the structure of the 
inner part of the disc. Here we describe this influence in the most crude approximation 
( a more accurate approach is discussed in the next subsection). The main effect of the 
secondary is to slow down the radial motion in the disc: v r <C tvo- To describe this effect 
quantitatively we divide the disc by two regions: 1) the outer region r > r in f ^> rj,, where the 
equations (44-47) are approximately valid, 2) the inner region ry, < r < rj„/, where the drift 
velocity (43) is small compared to the stationary value (46). These regions are separated by 
the 'influence' radius rj n //(t), and the dependence of rj n ^ on time is calculated below. When 
v r ~ 0, eq. (43) gives 

£(r,t) =C{t)r~^r. 48 

The functions C(t), and rj n //(t) are calculated using the laws of conservation of mass and 
angular momentum. The eq. of mass conservation has a form 

and the eq. of conservation of angular momentum has a form 

r r mfi(t) d o 3 o 

/ dr(-Xr 3 {l) = { V^£}_ . 50 

JO Ot A 

In both equations we formally extend the lower limit of integration to r = 0. In eq. (50) we 
neglect the small inflow of angular momentum from the outer (steady-state) region of the 
disc. After substitution of (47) into eqs. (49-50), we have 



1 1 2c- 1, M 
v = i 

infl I i i J a ' 

a + 1 Air 



2c+t 
r 2(a + : 

infl 4" a + 1 



3 r 2c + a 



CrB 1 ' = ~A^}kC a+1 . 52 
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The solution of eqs. (51-52) up to unimportant constant of integration has the form 

2c+a 

2c-l M 1 

C= (3fc)"^^(— ) {(3t)~c 53 

7T 

i M c 2c + a 2 a+i 
rw,= 3* - - (^T) ~ 54 

7T ZC — 1 



where 



c(2c + a) 2c + a N 



2(a + l) yv 2c-r 

Prior to the moment to the surface density in the disc is given by the steady-state solution 
(45), and after that moment, the surface density should be described by the solution (48, 
53). To reconcile these two solutions we set 

, _ W r o) 

where t„o( r o) = t u0 (r = r ). In this case eqs. (48), (53) can be written as 

S = S (r )(t/to)^(r/r )-^, 55 

and for the viscosity distribution we have 

V = ^(t/t )^(r/r y 1/2 , 56 
in 

where r\ = z/S. The expressions (55-56) are valid as long as the local viscous time t v = ^j- 
doesn't exceed the characteristic dynamical time td = 5 ~ t. We have 

t v /t ~ (r/r inf i)^, 57 

so our approximation breaks down at r ~ r in fi. Thus, to describe the disc structure at 
r ~ ri n fi we need a more accurate approach (see next subsection). 

Now we are able to calculate the evolution of the binary orbit. First, let us calculate the 
secondary radial velocity For that we use eq. (42), and then neglecting the second term 

dr b 67Tr b t](r b ,t) 



dt 

in the brackets we obtain ' 



dt m 



5s 



* * * In the case of a relatively massive secondary, the inner edge of the disc is determined by a few outer Lindbland resonances, 
and may be 2-3 times larger than the secondary orbital radius (e.g. Artymowicz & Lubow 1994). However we can use eq. (58) to 
calculate the radial velocity. This is connected with fact that the angular momentum outflow does not depend on the distance 
at all, see eq. (50). 
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and with help of eq. (56) eq. (57) can be rewritten in dimensionless form 

dy b S 



dr 2(3 ' 



where r = t/t , y b = (r(,/r ) 1//2 , and dimensionless parameter 



S=*5<1, 
m 



here M do = Mt u o(r ). The integration of eq. (59) gives 



5c+b 



2 



r 6 = (1 - 7 ^(r— - 1)) r , 60 

where 7 = p^ +b ^ ~ 1- The characteristic time-scale of the evolution of the secondary orbit 
t ev can be found from the condition r b (t ev ) = 0: 

tev = (r/S)~i&*t . 61 

It is instructive to compare the expression (60) with the Syer and Clarke result: 

n~ (1- (3 sc S^T)~r , 62 

1 + a 

where the dimensionless coefficient (3sc = 2 2+a p^+a) ~ 1- ^ can seen that for the typical 
values of the parameters a, b the evolution equation (62) underestimates the evolution time 
in comparison with equation (60). 

The dimensionless velocity (59) taken at the time t ev can be estimated as: 

1 dyb I n 4c 

_p ^ s^+z < 1. 
ar 

This inequality implies that we can treat the change of binary orbit as a slow effect with 
respect to the change of the disc configuration when t < t ev . On the other hand, let us 
look at the ratio of the disc mass contained inside the region influenced by the secondary 
(ri, < r < Tinfi) to the mass of the secondary. It can be easily seen that this ratio is small at 
the time £ ~ t ev : 

rn n fi(tev) 2(a+i) 
M d = 2tt / drrE(t ev ,r) ~ S^+^m < m, 63 

Jr b (t ev ) 

and this means we can use eq. (58) as a boundary condition to our problem. It is useful to 
express the evolutionary time t ev in terms of the accretion time t acc = ||. From eq. (61) we 
obtain 

tev = {lS)^^t acc < t acc . 64 

Say, for a = 2/3, b = 1, we have 

tev = (7*5*) ^ tacc} 65 
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and for a = 3/7, b = 15/14, we have 

t ev = {jS) 2/7 t acc . 66 

Finally let us estimate the luminosity of the disc. From the law of energy conservation it 
follows that the total luminosity of the disc L must be equal to the change of the secondary 
energy per unit of time: 

1 GmMdn _ GMM A i/ 2 f 

2 rl dt ~ r b (t) [ r b (ty [T) ' 

where is given by eq. (60). When the separation distance rj, is much larger than the 

gravitational radius of the primary, the "initial" value of luminosity (taken at the moment 
t = to) is much smaller than the "standard" value of luminosity of the steady state disc L st 

MM • erg 



'si 



10^M*M« 



r g s 



However, the luminosity grows with time due to evolution of inner part of the disc and the 
secondary orbit. 

4.2 Self-similar solutions 

Now we consider in detail the formal case of infinitely small binary separation distance 
(r& ~ 0). As we will see, in this case simple self-similar solutions of eq. (41) can be found 
f_ _ *| . These solutions are close to the exact solutions for the disc structure outside the binary 
orbit, and they allow us to find the disc structure at the intermediate scales r ~ r in fi. 

When considering the case r& = 0, we must specify the boundary conditions at r = 
and at infinity. We assume that there is no radial velocity at r = 0: v r (r = 0) = 0, and the 
disc structure tends to the steady solution with increase of r. We look for a solution of eq. 
(41) in the form 

a = t~ m F(0, 70 
f = D- X rt~ n , 71 

where we use the surface number density a = S/m p (m p is the proton mass), and the 
constant D is chosen to make the similarity variable (71) dimensionless. The powers n, m 

* * * Our approach is similar to the Lyubarskii and Shakura 1987 study of non-stationary disc accretion, and to Pringle's 1991 
calculation of evolution of external circumbinary disc of a given mass. 
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/2, 



ua), 72 



are specified by the condition that eq. (41), as well as the accretion rate 

£L = JL_ = 3r i/2 A f r v 

2/T 2ixm p dr ^ 

are expressed in terms of the variable £ only. Substituting eqs. (70) and (71) into eqs. (41) 

and (72), we obtain 



l + a b 
n = , m — -. 73 

c c 



Equation (41) takes the form 



with 



(3k)D b ~^e /2 ^e /2+b F 1+a + ZinS-fjrF + mF) = 0, 74 

a 

~ i 777 c 

£>=(3Jfe)*(— ) , 75 

where = m®k. The solution of eq. (74) should be separately analyzed for the general 
situation a ^ 0, and for the degenerate case a = 0. 

In the case a = we have (3k)D b ~ 2 = 1, Eq. (74) is linear, and can be easily solved. The 
general solution is 

,,£-1/2 a; -2(2-b) 

F = r b (^i + C 2 / dxe 15=S)3-). 76 



o 

The first term in the brackets represents the stationary solution (45) provided C\ — \. We 
have 

F « C 2 r b " 1/2 , 77 
at £ — > 0, and the comparison of eqs. (70,71), (77), and (48), (53) gives C 2 = (3. 

In the case a ^ 0, we can eliminate dimensional factor in the eq. (74) by redefining of 
the function F 

F = (3k)~«D^F. 78 

Substituting (73) (77) into (74), we obtain: 

d, 1/2 d^ 1/2+bpi+a + i {b p +{1 + a) ^d p = Q 7g 

at; a4 c a£ 

In the limit £ 3> 1 the solution of eq. (79) is close to the stationary solution 

F = £~^, 80 

and in the opposite limit ( C 1 we have 

b+l/2 

F = CC^ r - 81 
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Figure 4. Wc plot the numerical solution of eq. (79) F(£) as a function of the similarity variable £ (solid curves). Curve 1 
corresponds to the case a = 2/3, 6 = 1 (Thompson opacity dominates), and the curve 2 corresponds to the case a = 3/7, 
b = 15/14 (the opacity is determined by the free-free processes). The dashed lines represent the analytical asymptotic of that 
solutions at f — ► (see eq. (81)), and the dotted lines represent the stationary solutions of eq. (79) (see eq. (80)) 



It is evident from eqs. (70, 71, 81) that the asymptotic (81) corresponds to the approximate 
solution (48, 53) if we set again C = [3. Thus eqs. (45) and (48, 53) give asymptotic limits to 
the self-similar solution of eq. (41). The intermediate scale £ ~ 1 separates the asymptotic 
(80,81), and the corresponding radius r*(£ = l,t) is of order of r in fi. 

The numerical solution to eq. (79) is presented in Fig. 4. One can see from this Fig., 
that the solution of this equation approximated with high accuracy by the expression (81) 
at £ < 1, and by the expression (80) at £ > 1. We will use this fact below comparing our 
analytical estimates with the results of the numerical computations. 

4.3 Numerical model 

We solve numerically equation (41), and compare the results of our calculations with ana- 
lytical estimates. For numerical work we rewrite eq. (41) in terms of new radial coordinate 
y = (r/ro) 1 ^ 2 , dimensionless time r = t/t , and dimensionless surface density £ = S/S (r ): 




To solve this equation, we must specify the initial and boundary conditions. We assume that 
the disc surface density is equal to its steady state value (45) at some point y out ^> 1, and we 
take y ou t = 10 2 in our calculations. Since we are interested in the effect of evolution of the 
disc and the secondary orbit, the detailed description of the angular momentum exchange 





82 
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between the secondary and the disc is unimportant for us, and to obtain a boundary condition 
at y — yb we substitute the expression for the disc gas radial velocity (43) in eq. (42) and 
use the resulting equation as the boundary condition to our problem. The eq. (43) allows 
us to calculate the change of the secondary orbit. The initial surface density distribution 
£(£ = 0, y) is determined by eq. (45) in the range 2 < y < y out , and we use the appropriate 
polynomial function in the range 1 < y < 2 to match the steady state distribution and 
the inner boundary condition at the moment t — 0. The detailed form of this function is 

1 /2 

unimportant for us. The position of the secondary yi, = (r&/r ) is changing with time, 
and in our numerical scheme it is convenient to use another radial coordinate y which is 
stretching with change of the inner boundary with respect to coordinate r: 

~ (yout - l)y + y ou t(l -yb(t)) SQ 

y 77\ , so 

y ou t - yb(t) 

and y(yb(t)) = 1, y(y ou t) = yout- We use uniform mesh in coordinates (r, y) with 10 4 grid 
elements in y-direction, and employ a standard implicit scheme. To check the relaxation of 
numerical solution to the self-similar asymptote, we integrate eq. (82) with respect to time 
up to the time t out = 2t ev . When the orbital distance approaches the value yi n = 10~ 2 we 
"turn off" the radial evolution of the secondary. 

The results of numerical calculations are presented at Figs. (5-7). In Figs. 5,6 we show 
the calculated density profile taken in different moments of time, and for different param- 
eters of our model. One can see from these Figs, that the surface density of the disc is 
well approximated by the analytical expressions (45), (55). In Fig. 7 we show the depen- 
dence of the secondary orbital distance on time. The analytical curves give again very good 
approximation to the numerical results. 

5 ASTROPHYSICAL CONSEQUENCES AND CONCLUSIONS 

Now let us summarize the main results of our paper. We considered the interaction between 
a massive binary system and an accretion disc, paying special attention to the case when 
the secondary mass is larger than the disc mass inside the binary orbit. We showed that 
the joint evolution of the binary and disc could be divided into several successive stages. 
When the binary orbital plane is inclined with respect to plane of the disc, the gravitational 
field of the binary causes non-axisymmetric "twisted" distortions of an initially planar disc, 
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Figure 5. The dimcnsionlcss surface density £/Eo(ro) is plotted against dimensionless radius r/ro- The parameters a = 2/3 
and 6=1 correspond to the case of Thompson opacity "domination" . The different plots correspond to the different values of 
the parameter S: a) top left plot -S = 10 — , b) top right plot-S = 10 — 2 , c) bottom left plot -S = 10 — 3 , d) bottom right plot- 
s' = 10" 4 . The solid curves, in order of increasing of inner boundary, are made at times t = 10 2 t ev , t = 10 1 t e u and t = 2t ev . 
The dotted curves are determined by eqs. (45,55). In asymptotics the solutions have power-law form: E(r — » oo) ~ r -3 ' 5 , and 
E(r -> 0) ~ r" 9 / 10 . 

and the inner part of the disc lies in the orbital plane as a result of the evolution of these 
distortions. Since the orientation of outer part of the disc is controlled by the disc formation 
processes, after the alignment of the inner part of the disc a quasi- stationary twisted disc is 
formed. We calculate the radial dependence of the inclination angle of the twisted disc, and 
find the dependence of the characteristic alignment scale on viscosity and mass ratio of the 
binary components. It is shown that in the low-viscosity limit the alignment scale does not 
depend on the value of the viscosity parameter a. 

Assuming that the inner part of the disc lies in the binary orbital plane, rotates in the 
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Figure 6. The same as the plot b) of Fig. 5, but now a = 3/7, b = 15/14 (the opacity is determined by free-free processes). 
The general behavior of the curves is similar to the case b), but the asymptotical slope is slightly steeper: E(r — > oo) ~ r -3 / 4 , 
and E(r -> 0) ~ r" 11 / 10 
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Figure 7. The dimensionless binary separation distance n/ro is plotted against the time t/to. The parameters of the left plot 
are the same as in the plot b) of Fig. 5, and the parameters of the right plot are the same as in Fig. 6. The dashed curves are 
analytical expressions given by eq. (60), and the dotted curves represent Syer and Clarke result (62) 



same sense as the binary, we find the radial evolution of the disc structure and the evolution 
of the binary separation distance. In such a situation the disc gas cannot drift through the 
binary orbit, and is accumulated in some region outside the binary orbit. The structure of 
the disc in this region is close to the self-similar solution of the diffusion equation describing 
the viscous evolution of the disc surface density. The back reaction of the circumbinary disc 
pushes the secondary inwards towards the primary. We calculate analytically the dependence 
of the binary separation distance on time, and compare with results of numerical calculations. 
We find that the characteristic time scale of this evolution is determined by the ratio S of 
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the disc mass inside the initial binary orbit and the secondary mass m. The evolution time 
scale is always smaller than the "accretion" time scale t acc = |£. 

In the case of a supermassive black hole binary, evolution of the binary separation dis- 
tance is also caused by the emission of gravitational waves. To estimate the relative signif- 
icance of this effect, let us compare the total disc luminosity with orbital energy loss E gw 
due to emission of gravitational waves (Peters & Mathews 1963) f~ _ 

E gw = ^L*g 2 (^) 5 , 84 
5 r b 

where the "gravitational emission factor" = ^ ~ 3.6 x 10 59 erg/s. Comparing eqs. (69), 
and (84), we find that the interaction with the accretion disc is more important at scales 
larger than the characteristic "gravitational" scale r gw : 

r>r gw ^6- 10 2 ( T ^) 1/2 (M 8 M,)" V4 r 9 . 85 

As we have already discussed, the evolution of a supermassive black hole binary due to 
dynamical friction with the stars of the central stellar cluster can be unimportant at scales 
< Ipc, and therefore the secondary-disc interaction may provide the only possibility to pass 
through the "intermediate" scales r gw < r < Ipc. This fact may be critically important for 
formation of powerful sources of gravitational radiation due to coalescence of supermassive 
black holes. 

One uncertainty in our calculation seems to be connected with eccentricity evolution. In 
the course of binary evolution due to dynamical friction with the stars, the eccentricity seems 
not to grow significantly (see, e. g. Polnarev & Rees 1994 for analytic work and Quinlan Sz 
Hernquist 1997 for N-body simulations). However, the binary-disc interaction may induce 
some eccentricity (e.g. Artymowicz, 1992, Lin & Papaloizou, 1993). As we have mentioned, 
gravitational wave emission is much more effective for eccentric binaries, and a significant 
eccentricity can shorten the evolutionary time scale. 

* * * In the expression (84) we set the orbital eccentricity e = 0. For e ~ 1 the emission of gravitational waves is much more 
effective (Peters & Mathews 1963). 
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